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Abstract 


High  precision  Cold  Atom  Interferometers  (CAI)  are  in  development  to  supplement 
or  replace  conventional,  navigation  quality  inertial  measurement  units.  A  major  drawback 
of  the  atomic  interferometers  is  their  low  duty  cycle  and  sampling  rate,  caused  by  delays 
required  for  cooling  the  atoms  and  collecting  acceleration  and  angular  rate  measurements. 
A  method  is  herein  developed  for  inertial  navigation  by  integrating  highly  accurate,  low 
duty  cycle  CAI  measurements  with  high  bandwidth,  conventional  Inertial  Navigation 
System  (INS)  measurements.  A  fixed-lag  smoothing  algorithm  is  used  to  estimate  optimal 
acceleration  and  angular  rate  measurements  from  the  CAI  and  INS  data.  Given  current 
CAI  limitations,  simulation  results  demonstrate  nearly  50  percent  error  reduction  for  the 
enhanced  INS  compared  to  a  conventional,  unaided  INS.  When  the  conventional  INS 
position  error  was  increased  by  500  (m/hr),  the  50  percent  error  reduction  from  aiding 
was  maintained.  Increasing  the  conventional  INS  data  rate  fifteen-fold  while  maintaining 
a  1  Hz  CAI  sample  rate  leads  to  an  approximately  6  percent  increase  in  navigation 
error,  suggesting  that  the  CAI-aiding  algorithm  effectivity  is  only  slightly  influenced 
by  the  conventional  INS  data  rates.  A  five-fold  increase  of  the  CAI  measurement  rate 
shows  approximately  80  percent  reduction  in  navigation  error,  supporting  the  potential 
for  significant  performance  gains  in  the  near  future  from  advancements  in  cold  atom 
technology. 
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SIGNAL  PROCESSING  IN  COLD  ATOM  INTERFEROMETRY-BASED  INS 


I.  Introduction 


1.1  Research  Motivation 

Precision  navigation  is  critical  for  numerous  applications  in  both  civilian  and  defense 
sectors,  from  vehicle  guidance  to  accurate  target  location.  Historically,  passive  inertial 
navigation  was  the  dominant  guidance  system  used  for  naval  and  aviation  purposes,  but  lost 
favor  after  the  widespread  introduction  of  the  Global  Positioning  System  (GPS).  Recent 
advances  in  GPS  technology  now  provide  centimeter-level  measurements  for  users  in 
certain  applications.  Users  have  developed  a  dependency  upon  GPS,  a  weakness  which 
our  enemies  exploit  by  interfering  with  its  signals,  often  via  low  cost  methods.  In  addition 
to  deliberate  threats,  GPS  has  a  number  of  other  limitations.  Sub-meter  GPS  accuracy 
is  only  attainable  when  users’  receivers  have  uninterrupted  tracking  of  at  least  four  GPS 
satellites  [1].  Large  buildings,  canyons  and  indoor  environments  still  pose  problems  for 
users  attempting  to  acquire  and  track  GPS  signals.  In  light  of  these  limitations  and  the 
recognized  dependency,  interest  has  reverted  toward  alternative,  passive  sensors  for  aiding 
Inertial  Navigation  Systems  (INSs),  particularly  for  military  operations.  This  research 
will  investigate  the  feasibility  of  using  a  Cold  Atom  Interferometer  (CAI)  as  a  passive 
navigation  aid. 

Demonstrations  in  cold  atom  interferometry  technology  show  the  potential  for 
collecting  high  precision  measurements  in  the  fields  of  navigation,  gradiometry,  geophysics 
and  atomic  clocks  [2].  Numerous  universities  are  researching  hardware  designs 
and  improvements  for  cold  atom  interferometers,  including  Stanford  and  Ben-Gurian 
University  in  Israel  [3],  [4].  Stanford  University  has  demonstrated  the  performance  of  a 
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large,  1  (m^)  prototype  CAI  used  for  ground  navigation,  as  well  as  a  eompaet  CAI  for 
eolleeting  aerial  gradiometry  measurements.  The  CAI  performanee  eharaeteristies  used  in 
this  thesis  are  refereneed  from  the  Stanford  University  demonstrations  [5].  Ben-Gurian 
University  has  reeently  developed  a  CAI  fabrieated  ehip,  illustrated  in  Figure  1.1,  whieh 
opens  the  possibilities  for  eompaet,  affordable  applieations  in  the  near  future  [4]. 


Figure  1.1:  Atom  Chip  [6] 


1.1.1  Inertial  Navigation. 

Inertial  navigation  is  the  proeess  of  measuring  position,  veloeity  and  orientation 
through  use  of  dead-reekoning  inertial  sensors.  Speeifieally,  inertial  navigation  systems 
use  aeeeleration  and  angular  rotation  rate  measurements  eolleeted  from  aeeelerometers 
and  gyroseopes  (gyros),  and  then  integrate  these  measurements  to  determine  position, 
veloeity  and  attitude  [7] .  The  measured  aeeeleration  is  the  sum  of  linear  aeeeleration  and 


2 


gravity  along  the  accelerometer’s  input  axis  of  measurement.  Integration  of  the  gyroscopes’ 
rotation  rates  provides  angular  orientation  so  that  the  acceleration  vector  may  be  resolved 
in  a  reference  frame.  Typical  strapdown  inertial  navigation  systems  are  thus  composed 
of  three  accelerometers  and  three  gyroscopes  to  collect  measurements  from  the  three 
orthogonal  axes.  The  strapdown  INS  is  a  digital  alternative  to  the  mechanical  gimbal  and 
is  used  predominately  in  aerospace  navigation  because  of  its  improved  reliability,  lighter 
weight,  and  reduced  costs.  Ring  laser  or  fiber  optic  interferometers  act  as  rate  gyroscopes 
and  accelerometers  in  the  strap  down  systems  and  can  collect  measurements  on  the  order  of 
300  Hz.  The  INS  is  also  a  passive  system,  which  makes  its  operation  impervious  to  external 
denial  techniques  and  makes  them  very  useful  for  military  applications  against  jamming 
or  spoofing  threats  [7].  Despite  the  passive  nature  and  high  sampling  rates,  small  sensor 
errors,  which  are  integrated  as  part  of  the  inertial  measurements,  cause  the  INS  accuracy 
to  degrade  and  become  unreliable  over  time.  Thus,  the  INS  measurements  require  aiding 
from  other  sensors,  typically  GPS,  to  bound  error  growth  [7]. 

1.1.2  Cold  Atom  Interferometry. 

Conceptually,  the  cold  atom  interferometer  design  is  similar  to  an  optical  interferom¬ 
eter  used  in  both  ring  laser  and  fiber  optic  gyros.  Using  reflectors  and  beam  splitters,  a 
ring  laser  interferometer  guides  counter-propagating  light  beams  and  then  collects  rotation 
information  from  the  resultant  interference  patterns  as  illustrated  in  Figure  1.2.  Similarly, 
atomic  ’matter’  waves  are  split  using  either  material  structures  or  light  fields  [8]  as  shown 
by  the  n  and  |  beams  shown  in  Figure  1.2.  The  atom  waves  are  then  guided  using  ei¬ 
ther  lasers  or  magnetic  fields  [4].  The  wavelength  of  the  atomic  wave.  A,  is  related  to  the 
momentum  p  of  the  atom  according  to  the  quantum  mechanics  relationship  as  follows  [8]: 


p  mv 


where  h  is  Plank’s  constant,  m  is  the  atom’s  mass  and  v  is  its  velocity.  In  general,  CAI 
physics  is  based  upon  de  Broglie’s  wave-particle  duality  hypothesis  which  states  that  at 
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the  quantum  level,  matter  exhibits  wave-like  properties  [5].  In  the  CAI,  cooled  atoms  are 
guided  through  a  vacuum  and  their  resultant  acceleration-induced  phase  shifts  are  measured 
via  a  diffraction  grating.  The  phase  shift  lS.(piight  is  analogous  to  the  Sagnac  effect  in  the 
optical  interferometers  and  can  be  described  by  [5] 

A7TQ.h 

^^light  ~  (1*2) 

where  Q  is  the  angular  rate  of  the  instrument,  c  is  the  speed  of  light,  and  A  is  the  enclosed 
area  of  the  interferometer’s  optical  path.  The  phase  shifts  cause  interference  patterns  on 
the  diffraction  grating  from  which  phase  information  is  collected  and  used  to  determine 
angular  rotation.  Additionally,  because  of  their  particle  nature,  the  atoms  may  be  treated  as 
inertial  masses  and  their  movement  is  used  to  determine  the  CAI  system’s  acceleration. 

Despite  functional  similarities,  the  CATs  atom  wave  interferometer  offers  the 
advantage  of  higher  sensitivity  to  inertial  forces  as  compared  to  the  conventional,  optical 
interferometer  in  strapdown  INS.  The  group  velocity  of  the  atoms  is  relatively  small 
compared  to  the  speed  of  light,  so  very  small  accelerations  will  still  induce  a  significant 
change  in  the  atom  beam’s  path,  causing  measureable  phase  shifts  in  the  interference 
signal.  If  A  is  held  constant,  theoretically  CAI  is  more  sensitive  to  rotation  than  an  optical 
interferometer  by  a  factor  of  10  since  atomic  mass  is  greater  than  the  light  photon’s  relative 
mass  [3].  This  is  expressed  in  the  following  equation,  where  6(f>ught  was  defined  in  (1.2), 
and  6(patom  is  the  atom  phase  shift: 

^4^atom  _  Vn  ^  10^^ 

Sought  h/(^Ac^ 

Evaporative  cooling  techniques  allow  for  precise  control  of  the  velocity  and  position 
of  the  atoms  in  the  atomic  wave.  The  low  temperatures  necessary  to  cool  the  atoms  also 
improve  the  CATs  signal-to-noise  ratio  and  its  sensitivity  to  external  fields  [8]. 

However,  there  are  some  factors  which  reduce  the  CAI  sensitivity,  including 
bandwidth  and  duty  cycle.  Due  to  low  bandwidth  associated  with  high  precision  sensors. 
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Optical  Interferometer 


Figure  1.2:  Comparison  of  Optical  and  Atom  Interferometer  (Mach-Zehnder  type 
configuration). 


large  dynamics  may  cause  the  path  of  the  atom  particles  to  exceed  the  sensor’s  detection 
threshold.  Thus,  a  high  dynamic  environment  may  produce  unreliable  measurements  or 
CAI  system  failure.  Additionally,  the  time  duration  for  collecting  a  CAI  measurement  is 
on  the  order  of  1  Hz  [8],  which  leads  to  a  loss  in  collected  information  as  compared  to  the 
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typical  50-300  Hz  sample  availability  provided  by  conventional  strap-down  INS.  The  main 
themes  for  this  research  are  both  understanding  and  compensating  for  this  latter  issue. 

1.2  Problem  Definition 

The  objective  of  this  thesis  is  to  explore  Inertial  Navigation  using  CAI-based 
acceleration  and  angular  rate  measurements.  While  CAI-based  measurements  are  very 
accurate,  the  rate  at  which  they  can  be  obtained  is  relatively  low,  e.g.  1  Hz.  The  low 
duty  cycle  is  caused  by  delays  required  for  cooling  the  atoms  and  collecting  acceleration 
and  angular  rate  measurements  [5],  [8].  Less  accurate,  conventional  INS  measurements, 
which  are  available  at  a  much  higher  rate,  might  be  used  to  fill  the  temporal  gap  in 
CAI  measurements.  In  addition,  the  CAI  measurements  represent  an  average  of  the  true 
acceleration  and  angular  rate  during  the  collection  phase  of  the  duty  cycle.  Therefore,  the 
CAI  measurement  may  be  combined  with  the  conventional  INS  data  during  this  phase 
to  reduce  the  overall  navigation  error.  For  a  CAI  sensor  with  50  percent  duty  cycle, 
the  conventional  INS  measurements  would  be  used  for  navigation  during  the  first  half 
of  the  duty  cycle  where  CAI  measurements  are  unavailable.  Subsequently,  a  weighted 
combination  of  1)  the  conventional  INS-provided  multiple  measurements,  and  2)  the  single, 
CAI-based  acceleration  and  angular  rate  measurements  would  be  used  for  the  second  half 
of  the  duty  cycle. 

1.3  Related  Research 

1.3.1  Multiple  INS  Integration. 

Omr,  Georgy  and  Noreldin  propose  a  method  for  integrating  the  accelerometer 
measurements  from  multiple  MEM  sensors  for  pedestrians  in  [9].  The  authors  show  that 
Micro  Electro-Mechanical  System  (MEMS)  measurements  from  multiple  devices  carried 
on  different  parts  of  the  user  may  be  collected  and  combined  wirelessly.  Associated 
covariances  for  the  position  and  velocity  measurements  are  used  to  calculate  the  best 
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estimate  for  a  unified  navigation  solution,  or  conversely,  used  to  update  all  the  individual 
MEM’s  measurement  solutions.  The  covariance  weighting  algorithm  requires  less 
computing  power  than  a  centralized  Kalman  filter  because  it  only  updates  the  solutions.  In 
other  words,  propagation  of  an  internal  Kalman  filter  model  is  unnecessary.  This  method 
allows  for  integrating  an  unrestricted  number  of  inertial  navigation  sensors,  which  the 
authors  note  is  a  challenge  in  centralized  Kalman  filters. 

This  thesis  differs  from  Omr’s  research  in  that  this  work  considers  more  complex 
inertial  navigation  sensors,  corresponding  measurements,  and  purpose.  The  cold  atom 
interferometer  and  INS  sensor  are  much  more  accurate  than  MEMs  and  provide  angular 
rate  measurements  in  addition  to  acceleration  measurements.  Thus,  while  Omr’s  research  is 
limited  to  aiding  position  and  velocity,  this  thesis  also  proposes  methods  for  aiding  attitude. 
Additionally,  this  thesis  accounts  for  the  problems  posed  by  integrating  measurements  from 
sources  which  operate  at  different  sample  rates  and  duty  cycles.  Einally,  while  Omr  et 
al  focus  on  pedestrian  navigation,  this  thesis  applies  multiple  sensor  integration  to  three 
dimensional  aircraft  navigation. 

1.3.2  CAI  Integration  with  Conventional  INS. 

Jekeli  [5]  proposes  a  model  for  the  cold  atom  interferometer  which,  with  some 
reasonable  assumptions,  can  be  viewed  as  analogous  to  the  conventional  INS.  Using  bias 
and  drift  parameters  associated  with  the  accelerometers  and  gyroscope,  he  demonstrates 
the  CAI  inertial  sensor’s  theoretical  performance.  The  acceleration  and  gyroscope 
measurement  models  presented  in  Jekeli’s  research  [5],  and  substantiated  in  [10],  are  used 
for  the  CAI  and  conventional  INS  measurements  in  this  thesis. 

Canciani  and  Raquet  propose  three  filter  frameworks  to  integrate  a  CAI  with  a 
conventional,  navigation-grade  INS  [11].  The  first  framework  corrects  the  INS  acceleration 
and  attitude  measurements  and  when  available,  mechanizes  off  the  CAI  measurements. 
The  second  framework  always  uses  the  conventional  INS  in  the  mechanization  equations. 
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but  uses  the  CAI  delta-velocities  and  delta-thetas  to  correct  the  navigation  solution  at  the 
position  level.  The  third  framework  proposes  integrating  CAI  with  GPS  measurements 
in  order  to  compensate  for  GPS  outages.  When  outage  times  were  very  short,  the  first 
framework  showed  the  best  results,  while  the  second  framework  showed  the  best  results 
during  the  longer,  high-g  CAI  outages.  The  simulation  tests  assumed  that  the  CAI 
measurements  would  be  available  consistently  until  outages  occurred;  for  example,  a  duty 
cycle  of  80  percent  was  modeled  as  a  square  wave  with  periodic  2  minute  outages. 

While  Canciani’s  research  focused  on  the  CAI  dropouts  and  outages,  this  thesis 
investigates  the  effects  from  the  CATs  low  duty  cycle,  sampling  rates,  and  the  analysis 
methods  to  compensate  for  them.  This  thesis  makes  measurement  corrections  on  the 
acceleration  and  angular  rate  level,  similar  to  the  first  framework.  However,  Kalman  filter 
modeling  is  not  used;  rather,  this  thesis  implements  an  estimation  algorithm  similar  to  a 
fixed-lag  smoother. 

1.4  Scope  and  Assumptions 

This  research  uses  a  proof  of  concept  approach  to  model  the  CAI  and  INS  errors  and 
performance  characteristics.  The  acceleration  and  gyroscope  error  models  are  composed 
of  a  bias  and  discrete  measurement  drift  based  upon  a  continuous-time,  stochastic  process. 
The  analytical  development  of  the  latter  is  made  in  Appendix  A.  A  linear  relationship  is 
assumed  between  the  sensor’s  drift  induced  measurement  variances  and  bias  variances, 
which  allows  for  calibration  using  the  Lyapunov  covariance  equation  in  Appendix  B. 
Furthermore,  the  Earth  model  is  simplified  to  a  flat,  non-rotational  environment  to  simplify 
the  mechanization  equations.  This  assumption  neglects  the  effects  of  Coriolis  acceleration, 
Schuler  cycle,  transport  rate  and  other  frame  translational  errors  introduced  from  the 
Earth’s  curvature  and  rotation.  Altitude  aiding  from  barometric  data  or  radar  altimeter 
measurements  is  neglected,  although  it  would  be  generally  present  in  conventional, 
navigation-grade  INS. 
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As  noted  earlier,  this  research  focuses  on  the  limitations  of  CAI  aiding  pertaining  to  the 
low  duty  cycle  and  data  rate.  Additional  effects  from  high  dynamic  dropouts  are  neglected, 
although  it  is  assumed  that  an  integrity  monitor  would  exist  in  a  fielded  implementation  to 
disregard  CAI  data  when  dynamics  exceeded  an  experimentally  determined  g-level. 

1.5  Methodology 

In  order  to  integrate  the  CAI  with  the  INS,  measurement  models  were  developed 
and  error  parameters  were  calibrated  to  meet  their  respective  theoretical  performances. 
These  error  parameters  were  used  to  create  stochastic  random  variables  which  corrupted  a 
simulated  truth.  The  availability  of  the  CAI  and  conventional  INS  measurement  samples 
was  modeled.  An  analytical  method  was  developed  based  upon  this  model  to  combine 
the  measurement  samples  and  estimate  their  bias  and  drift  variances  for  each  periodic 
CAI  duty  cycle.  MATLAB  software  was  used  to  simulate  noisy  conventional  INS  and 
CAI  measurements  and  then  integrate  them  according  to  the  proposed  algorithm  methods. 
The  resultant  error  means  and  standard  deviations  were  tabulated  a  part  of  a  Monte  Carlo 
analysis. 

Chapter  2  presents  the  mathematical  background  necessary  to  understand  basic 
inertial  navigation,  error  modeling,  and  basic  Kalman  Filtering.  Chapter  3  develops  the 
aforementioned  measurement  models  and  the  proposed  algorithms  for  integrating  them. 
Chapter  4  presents  the  methods  used  to  statistically  analyze  the  simulation  results.  Error 
characterizations  were  made  by  evaluating  the  statistical  results  of  multiple,  hour-long 
runs  to  nullify  the  effects  of  potential  anomalous  behavior  in  single  runs.  Chapter  5 
presents  a  conclusion  of  the  data  analysis,  CAI  aiding  performance,  and  recommendations 
for  improvements.  Areas  of  future  research  are  also  suggested  based  upon  predicted 
technological  developments. 
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II.  Mathematical  Background 


This  chapter  presents  the  background  information  required  to  understand  the 
navigation  topics  and  mathematics  discussed  in  the  presented  research.  It  covers  basic 
notation,  reference  frames,  inertial  navigation  system  models  and  an  overview  of  Kalman 
filters. 

2.1  Notation 

Notation  used  consistently  is  described  below: 

Scalars:  denoted  as  upper  or  lower  case  italic  letters. 

Vectors:  denoted  as  lower  case  bold  font  letters  and  assumed  as  column  vectors  unless  oth¬ 
erwise  specified. 

Matrices:  denoted  as  upper  case  bold  font  letters  and  may  have  subscripts  specifying  rows 
and  columns. 

Transpose:  denoted  by  a  superscript  T. 

Estimated  variables:  denoted  by  the  hat  superscript  (e.g.  v). 

Reference  frame:  denoted  by  superscript  in  parentheses  (e.g.  is  a  scalar  in  the  ’a’ 
frame). 


2.2  Reference  Frames 

A  reference  frame  is  necessary  for  expressing  and  measuring  the  motion  of  bodies  in 
navigation.  Depending  upon  the  scope,  some  reference  frames  are  more  convenient  than 
others.  All  reference  frames  discussed  herein  use  three-dimensional  axes,  with  right  hand, 
orthonormal  basis. 


10 


2.2.1  True  Inertial  Frame. 


A  true  inertial  frame  is  a  theoretieal,  non-aeeelerating  frame  from  whieh  Newton  based 
his  first  law  of  motion.  Due  to  relativity,  this  frame  has  no  origin  and  is  stationary  relative 
to  absolute  time. 

2.2.2  Earth-Centered  Inertial  Frame  (i-frame). 

The  Earth-eentered  inertial  frame  is  defined  by  the  fixed  stars,  with  an  origin  at  the 
eenter  of  the  earth.  The  vertieal  axis  aligns  with  the  North  Pole  and  horizontal  x-axis 
and  y-axis  lie  along  the  equatorial  plane.  This  frame  is  non-rotating,  but  does  follow  the 
Earth’s  rotation  about  the  sun.  However,  beeause  the  Earth’s  orbit  about  the  sun  is  very 
slow  in  relation  to  navigation  about  the  earth,  the  i-frame  may  be  approximated  as  a  true 
inertial  frame. 

2.2.3  Earth-Centered  Earth-Fixed  Frame  (e-frame). 

As  its  name  implies,  the  Earth-Centered  Earth-Eixed  (ECEE)  frame  also  has  an  origin 
at  the  eenter  of  the  earth.  However,  its  horizontal  axes  move  to  eoineide  with  the  earth’s 
rotation.  The  x-axis  points  to  the  Greenwieh  Meridian  and  the  y-axis  points  to  90  degrees 
east  longitude  [7]. 

2.2. 4  Navigation  Frame  ( n-frame ). 

The  navigation  frame  origin  eoineides  with  the  loeation  of  the  navigation  system  of 
interest.  The  vertieal  axis  points  downward  in  the  direetion  of  the  loeal  vertieal-defined 
as  the  gravity  veetor  for  a  given  loeation  on  the  Earth.  The  x-axis  and  y-axis  point  North 
and  East,  respeetively,  giving  rise  to  the  familiar  North-East-Down  (NED)  eonvention.  The 
navigation  system’s  motion  eauses  the  navigation  frame  to  rotate  at  the  transport  rate  a>e„. 

2.2.5  Body  Frame  (b-frame). 

The  body  frame  origin  is  eo-loeated  with  the  n-frame,  but  rigidly  attaehed  to  the 
vehiele.  That  is,  the  x,  y,  and  vertieal  z  axes  point  in  the  direetion  of  the  aireraft  nose, 
right  wing  and  bottom  of  the  aireraft,  respeetively.  The  axis  set  rotates  with  roll,  piteh. 
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and  yaw  angles  as  the  aircraft  moves.  It  is  important  to  note  that  in  strapdown  INS,  the 
accelerometer  and  gyroscope  measurements  are  initially  resolved  in  the  b-frame  and  then 
converted  to  the  navigation  frame  using  the  Euler  angles  via  Direction  Cosine  Matrices 
(DCM). 


2.3  Inertial  Navigation 

2.3.1  Mechanization  Equations. 

The  following  equations  are  derived  from  [11]  and  [7].  For  this  research,  the  reference 
frame  for  the  systems  of  interests  will  be  the  local  geographical  level,  i.e.  navigation  frame. 
The  inertial  navigation  system  accelerometers  measure  the  specific  force  vector  f,  which  is 
the  difference  between  the  inertial  accelerations,  and  gravity 

jKO  _  j;(0  _  g(0 


Rearranging,  the  acceleration  in  the  naviation  frame  with  respect  to  the  earth  is: 

=  +  (2.1) 

where  g*^”^  is  gravity  in  the  navigation  frame,  is  specific  force  in  the  navigation  frame, 
cof^  is  the  rotational  rate  of  the  Earth,  (Jg^  is  the  transport  rate,  and  the  term  (2a>[”^  +  (J^n)  is 
the  Coriolis  acceleration  caused  from  navigating  in  a  non-inertial,  rotating  reference  frame 
[12].  Integrating  (2.1)  gives  the  three  dimensional  velocity  vector 

=  vf  (2.2) 


To  simplify  analysis,  this  research  assumed  a  nonrotating  and  flat  Earth,  which  ignores 
the  effects  of  the  Coriolis  acceleration,  so  that  the  n-frame  =  i-frame.  The  inertial  system 
accelerometers  only  provide  measurements  in  the  body  frame,  so  a  DCM,  C|f^  is  used  to 
resolve  the  force  vector  into  the  navigation  frame 


=  C 


(«)(-)(*) 
b  ^nb 
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In  our  simplified  inertial  frame,  the  DCM  beeomes  simply  the  identity  matrix  and 


is  assumed  equivalent  to 


2W 

Figure  2.1:  Visualization  of  Loeal  Level  Plane  with  assoeiated  body  frame  and  navigation 
frame. 


2.3.2  Strapdown  Error  States. 

Sinee  a  non-rotating,  flat  Earth  is  assumed,  flight  in  the  loeal  level  plane  may  be 
visualized  as  Figure  2.1  where  the  speeifie  foree  veetor  f  is  eomprised  of  the  axis  speeifie 
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components  fy.,  and  of  the  i-frame  [13] 


fx.  =  fx,cose  +  fz^sine 

fy.  =  0 

fz.  =  -fx,sine  +  f^,cose  (2.3) 

where  9  is  the  aircraft  pitch  angle  and  fx^  and  are  the  axis  specific  force  components  in 
the  b-frame.  Perturbing  (2.3)  gives  the  specific  force  error  equations 

Sfxi  =  i-fxtSine  +  f^,cos9)6e  +  dfx.cosO  +  df^^sinO 
=  fz,59  +  6fx,cose6f^,sine 
5fy,  =  0 

5fzi  =  -{fx.cosO  +  f^^sine)5e  -  dfx.sinO  +  df^^cosO 

=  -fx,6e  -  6fx,sine6f^^cose 


The  errors  in  acceleration,  velocity,  and  angular  rate  are  found  by  perturbing  equations 
(2.1)  and  (2.2): 

6xi  =  6vx, 

Sji  =  0 
6zi  =  dv^. 

5vx.  =  Sfx,  =  +  coses  fx,  +  sines 

Svy,  =  0 

dv^,  =  d/^,  =  -fxfe  -  sines  fx,  +  coseSf,, 
se  =  Sto  (2.4) 

where  di,,  dy,-,  and  Szi  are  the  axis  specific  velocity  errors;  dv;^;,  dv^,,,  and  dv^,  are  the 
axis  specific  acceleration  errors;  and  Se  is  the  angular  rate  error. 
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If  we  further  assume  wings  level,  constant  altitude  flight,  then  the  vertical  force 
simplifies  to  the  gravitational  acceleration,  g;  the  horizontal  force  6fxj  simplifies  to  the 
horizontal  acceleration,  a,  and  6  =  0.  The  equations  from  (2.4)  are  combined  to  form  a 
9x9  error  state  vector  6x  as  follows: 


6x  = 


dp  dv  dT 


where  dp,  dv  are  the  position  and  velocity  errors  vectors,  and  d'P  is  a  3  x  1  vector  of  angular 
errors  due  to  misalignment  from  the  navigation  frame,  that  is 


d'P  = 


cMn)  _  p(fo) 


=  -dC 


(n) 


-6i//  -69  -6(f) 


The  navigation  state  error  dx  equation  (with  driving  inputs  df^^^  and  is 


dx  =  Mdx  +  r 


df(^) 


6(1) 


where 


0 

I 

0 

0 

0 

0 

0 

p(n) 

,  r  = 

p(n) 

0 

0 

0 

0 

0 

f^(n) 

In  wings  level  flight,  is  the  skew  symmetric  form  of  the  specific  force  vector 


0  g  0 

-g  0  -a 
0  a  0 


2.4  Kalman  Filter 

The  Kalman  filter  is  a  linear,  recursive,  optimal  data  processing  algorithm  used 
to  estimate  the  unknown  states  of  a  stochastic  differential  equation.  Using  current,  noisy 
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measurements  of  some  components  of  the  state,  and  a  state  estimate  based  on  the  previous 
measurements,  the  Kalman  filter  calculates  the  current,  statistically  optimal  states.  The 
predicted  current  state  is  more  accurate  than  the  measured  state,  which  is  affected  by 
sensor  noise.  Unlike  batch  estimation  techniques,  storing  a  history  of  measurements 
or  estimates  is  unnecessary  because  the  statistical  information  is  contained  within  the 
Kalman  filter’s  recursive  state  estimates  and  associated  covariance.  The  Kalman  filtering 
algorithm  is  based  upon  the  following  key  assumptions:  the  linear  mathematical  model 
is  an  accurate  representation  of  the  true  system,  all  the  random  process  errors  in  the 
model  and  measurements  are  white  Gaussian  noises  and  the  measurement  observations 
are  uncorrelated  [14],  [15].  It  should  be  noted  that  no  true  process  is  perfectly  described  by 
a  linear  model,  thus  filter  tuning  and  modification  are  sometimes  necessary  to  achieve  the 
best  solution. 

2.4.1  State  Model. 

The  following  Kalman  filter  equations  closely  follow  [14].  The  linear  system  model 
satisfies  the  stochastic  differential  equation: 

x{t)  =  Ax(t)  +  Bu(t)  +  Gw(t)  (2.5) 

with  state  vector  x(t),  control  input  u(0  and  system  process  noise  w(0-  The  A  matrix  gives 
the  mathematical  coupling  between  states  and  the  G  matrix  relates  the  system  dynamics 
noise  to  the  states.  The  dynamics  noise  w(0  is  modeled  as  zero-mean  white  Gaussian  noise 
with  a  noise  strength  Q  described  by 


E[w(0]  =  0 

E[w(0w(t')^]  =  Q{t)6{t  -  t') 
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Integrating  equation  (2.5)  gives  the  discrete-time  form  solution  [14]: 


x,-+i 


=  0(At)X;  +  Wrf, 


where  is  a  Brownian  motion  process  with  dispersion  Q,  and  O  is  the  discrete-time 
state  transition  matrix  [14].  That  is, 

ti)  =  O(A0  = 


yidi  is  a  zero-mean,  uncorrelated  (i.e.  white  Gaussian)  process  with  discretized  noise 
strength  Q^.: 


Qd.  =  E[Wrf,wJ] 

0(t,+i ,  T)G(r)Q(T)G(T)'^0(t,+i , 

2.4.2  Measurement  Model. 

The  Kalman  filter  accepts  discrete-time  measurements/observations,  z,  from  external 
sensors  and  incorporates  them  as  part  of  its  optimal  estimate.  Recall  that  the  measurements 
must  be  uncorrelated  both  in  time  and  with  other  measurements.  The  measurements  must 
be  of  the  linear  form 

Z;  =  H;X;  -t  V; 


The  H  matrix  is  the  mathematical  coupling  betwen  the  system  states  in  x  with  the 
measurements  in  vector  z.  The  measurement  noise  vector,  v  represents  zero-mean  white 
Gaussian  noise  with  strength  R: 


E[v,-v}]  = 


R,  for  i  =  j 
0  for  i  4^  j 
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2.4.3  Kalman  Filter  Algorithm. 

As  discussed  above,  the  Kalman  filter  is  developed  using  a  model  of  the  dynamic 
system,  where  a  state  transition  matrix  O,,  input  matrix  B,,  and  noise  vector  w  with  strength 
Q,  are  used  to  determine  the  state  of  the  model.  Given  knowledge  of  the  initial  state  and 
covariance  estimates  of  the  system,  updates  to  state  estimates  can  be  calculated.  State 
and  covariance  updates  are  designated  by  x  and  P,  respectively.  The  algorithm  estimation 
is  comprised  of  two  phases:  propagation  and  updating.  In  the  propagation  phase,  the 
current  state  estimates  are  transitioned  forward  in  time  according  to  the  system  model. 
In  the  update  phase,  external  measurement  information  is  incorporated  and  the  Kalman 
filter  creates  a  new  optimal  estimate  for  the  states  using  this  information.  State  and 
covariance  estimates  updated  by  measurement  z„  are  designated  by  a  superscripted  plus 
sign,  and  P^,  respectively.  Similarly,  state  and  covariance  estimates  propagated  forward 
in  time  without  measurement  updates  are  designated  by  a  superscripted  minus  sign,  x“  and 
P“,  respectively.  Since  the  system  and  measurement  models  are  stochastic,  all  states  and 
measurements  are  treated  as  Gaussian  random  processes  [12].  One  of  the  unique  properties 
of  Gaussian  random  variables  is  that  they  can  be  completely  described  by  their  first  two 
moments.  Therefore,  only  the  states’  mean  and  covariance  are  calculated,  propagated  and 
updated  during  the  Kalman  filter  phases. 

The  Kalman  filter  algorithm  begins  with  initial  conditions  xjlj  and  P^j  as  inputs  to 
the  propagation  phase.  The  propagation  routine  steps  from  time  /  -  1  to  i.  0;_i,  Q,_i,  B,_i, 
and  control  input  u  are  used  to  calculate  x“  and  P^  according  to 

xj  =  +  B,_iu,_i 

P7  =  0,_iPti®L+G^QGa 

x“  and  P^  are  now  the  state  and  covariance  estimates  propagated  forward  to  the  next 
step.  If  no  updated  measurements  are  available,  then  the  current  state  is  reiterated  through 
the  measurement  routine.  If  a  measurement  z  is  available,  the  update  phase  is  initiated.  The 
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state  and  covariance  estimates  from  the  propagation  phase  are  used  as  inputs  to  the  update 
phase.  Given  the  modeled  measurement  matrix  H,  and  measurement  noise  strength  R,,  the 
optimal  Kalman  gain  K,  is  calculated  as  follows  [14]: 

The  Kalman  gain  is  then  used  to  determine  the  amount  of  influence  that  the  residual 
r,  discussed  in  detail  below,  will  have  on  the  updated  state  estimate  x^.  A  new  covariance 
estimate  Pf  is  also  calculated  using  the  optimal  Kalman  gain  according  to  [14]: 


xl  =  X;  +  K;r; 

(2.6) 

p;  =  P“  -  K,H,p: 

(2.7) 

r,  =  z,  -  H;x7 

(2.8) 

Both  xl  and  can  be  output  as  the  updated  estimates  and  utilized  in  the  state 
propagation  estimates.  Prior  to  entering  the  propagation  routine,  the  time  index  is  reset 
and  and  P^  become  xj^^  and  Pl_j.  A  summary  of  the  Kalman  filter  algorithm  is  shown 
in  Figure  2.2. 

In  a  modeling  and  simulation  setting,  a  truth  model  is  used  to  refine  the  Kalman  filter 
and  improve  its  outputs.  The  difference  between  the  model’s  truth  and  the  Kalman  filter  is 
known  as  the  error  [14]. 

el  =  Xi  -  xl  (2.9) 

Assuming  the  Kalman  filter  state  estimates  are  Gaussian,  then  the  error  vector  must  also  be 
Gaussian  and  completely  described  by  its  mean  and  covariance  [14],  [16]. 

E[ef|z;]  =x,+  -xf  =  0  (2.10) 

E[elef\z]=Pl  (2.11) 
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Reset  Time  Index 

(t  =  t  -  1) 


Figure  2.2:  Kalman  Filter  Process  [15] 


The  error  covariance  is  completely  independent  of  the  measurements,  but  are  dependent 
on  the  type  and  quality  of  the  measurement  [14].  For  analysis  in  this  thesis,  the  estimation 
error  can  be  used  to  quantify  the  estimation  algorithm’s  effectiveness. 

In  reality,  the  truth  data  is  not  known  and  must  be  measured.  Noise  is  often  introduced 
with  actual  measurements  and  must  be  weighed  against  the  estimate  of  the  next  state. 
A  residual  is  a  comparison  of  the  Kalman  filter  output  to  the  measurement  states.  The 
residual  is  defined  as  the  difference  between  the  measurement  and  the  state  estimate  prior 
to  measurement  update,  as  shown  in  equation  (2.8).  Residuals  associated  with  the  present 
state  are  independent  of  past  residuals  and  are  described  as  white,  Gaussian,  and  zero  mean 
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[16].  Like  error,  the  residual  can  also  be  described  by  its  mean  and  covariance  as  follows 
[14]. 


E[r,]=0  (2.12) 

E[r,rT]  =  H,P7HT  +  R,  (2.13) 

Given  a  nominally  operating  Kalman  filter,  residuals  can  be  statistically  monitored  for 
changes,  such  as  a  bias  shift.  Significant  changes  typically  indicate  a  malfunction  or 
anomaly  in  the  measurement,  for  example  a  failed  sensor  or  “noisy”  failure  [14].  Figure  2.3 
illustrates  residual  monitoring  and  sensor  failure  detection  due  to  residual  bias  shift  (top) 
and  residual  strength  increase  (bottom). 


BIAS  SHIFT  FAILURE 


“NOISY’  FAILURE 


I 

Figure  2.3:  Residual  Monitoring  and  Sensor  Failure  Detection  [14] 
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As  opposed  to  simulation  studies,  true  states  cannot  be  measured  directly  in  practice. 
Therefore,  it  is  difficult  to  determine  if  the  Kalman  filter  is  accurately  estimating  the  states 
and  whether  the  model  is  an  accurate  representation  of  reality.  This  uncertainty  is  mitigated 
by  rigorous  testing  and  simulation  of  the  model  by  generating  a  simulated  truth  based  on 
knowledge  of  how  the  system  behaves  dynamically  and  how  measurements  are  reported. 
Chapter  3  will  detail  specifics  of  these  system  dynamics. 


21 


III.  Algorithm  Development  and  Model  Creation 


This  chapter  develops  the  INS  and  CAI  accelerometer  and  gyroscope  sensor  models, 
the  algorithm  model  used  to  integrate  the  CAI  and  INS  measurements,  and  the  navigation 
mechanization  equations  used  for  system  simulation. 


3.1  Postulate 

In  principle,  the  CAI-based  acceleration  measurement  provided  at  the  discrete 
time  instant  k  at  the  completion  of  the  kth  duty  cycle  is  the  average  acceleration  during  the 
CAI-based  measurement  interval.  In  other  words,  it  is  the  area  under  the  true  acceleration 
divided  by  the  CAI  duty  cycle  length,  crAT : 


iCAI) 


Ja-a)A 


u{t)dt 


(3.1) 


(k-a)AT 


where  a  is  the  CAI  accelerometer  or  gyro  duty  cycle  fraction,  AT  is  the  sampling  time 
of  the  CAI-based  acceleration  measurements,  and  u(t)  is  the  true  acceleration.  Figure  3.1 
illustrates  the  derivation  for  (3.1). 

During  the  CAI  accelerometer  or  gyro  initiation  time  interval  [(A:  -  l)Ar,  {k  -  a) AT) 
in  duty  cycle  k,  there  are  N  -  n  conventional  INS  acceleration  measurements  Um,  available, 
I  =  0, 1, ...,  A-n- 1.  In  addition,  during  the  duty  cycle  interval  [(A:  -  Q')Ar,  MT),  where  the 
CAI  accelerometer  performs  an  acceleration  measurement,  there  are  n  +  \  measurements 
of  the  conventional  INS,  for  /  =  N  -  n,N  -  n  +  l,  ...N .  A  is  the  total  number  of  acceleration 
measurements  from  the  conventional  INS  during  the  CAI-based  accelerometer’s  sampling 
time  interval  AT .  Thus, 

(»1) 

ot 


where  5t  is  the  short  sampling  time  of  the  conventional  INS  accelerometer.  The  number 
of  conventional  INS  acceleration  measurements  taken  while  the  CAI  accelerometer  is  also 
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collecting  a  measurement  is  n  +  1,  as  shown  in  Figure  3.1.  Thus, 


A 

n  = 


aAT 


M(t) 


Figure  3.1:  Model  1:  Acceleration  over  one  duty  cycle  {k) 


The  sampling  intervals  AT  and  6t  of  the  CAI  and  conventional  accelerometers, 
respectively,  and  the  duty  cycle  fraction  a,  are  such  that  both  N  and  n  are  integers.  The 
area.  A,  under  the  acceleration  u(t)  during  the  duty  cycle  interval  [(k  -  a) AT  <t<  kAT)  is 
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approximated  using  the  trapezoidal  integration  rule,  as  shown  in  the  following  equation: 


^  ~  2  “  a)AT)  +  u((k  -  a)AT  +  5t)]  6t  +  ... 

+  -  [u(kAT  -  6t)  +  u{kAT)\  6t 


N-\ 

E 


+  >  M 

i=N-n-¥\  ) 


6t 


(3.2) 


where  m,-  =  u(i  ■  6t). 

The  CAI-based  accelerometer  and  gyro  measurements  are  subject  to  drift.  Each  of  the 
CAI  accelerometer-provided  acceleration  samples  has  a  small  error  caused  by  drift 
~  N(0,  The  quantification  of  the  standard  deviation  of  the  sampling  error  ctcai 
is  discussed  in  Appendix  A.  From  (3.1)  and  (3.2),  the  true  CAI  measurement  is  obtained: 


ACAI) 


6u 


(CAI) 


^kAT 

0:AT  J(k-a)AT 


u{t)dt 


N-l  \ 

+  ^  u 

i=N-n+l  J 


(3.3) 


Expanding  (3.3)  gives 


=  -u^-n  + 


m 


-l- ...  -I-  Mai_i  -I-  2^'^'  ^ 


(3.4) 


While  the  accleration  measurement  is  exclusively  discussed  herein,  the  same  also 
applies  to  the  measurement  of  the  aircraft’s  angular  rate  measured  by  the  CAI-based  gyros. 
In  particular,  ctcai  would  be  substituted  with  the  corresponding  drift  uncertainty  for  the  CAI- 
based  accelerometer  (cr),"'* j  and  the  CAI-based  gyroscope  discussed  in  Appendix  A. 


3.2  Signal  Processing 

Each  of  the  acceleration  and  angular  rate  samples  provided  by  the  conventional  INS 
accelerometers  and  gyros  have  a  constant  error  from  an  unknown,  random  bias  6uh  as 
well  as  error  from  drift  ~  A/'(0,  cr^^^);  obviously,  ctcai  «  o-,ns-  During  the  second 
half  of  each  duty  cycle  [N  -  n,  ...N],  there  are  n  -1-  1  acceleration  measurements  from  the 
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conventional  INS, 


u'Z-n  =  +  5uZl'' 


uZZ  =  u,  +  6ut  +  6uZ^^ 


(3.5) 


- ik) 

The  conventional  sensor’s  bias  6uh  is  recursively  estimated.  Thus,  the  estimate  6uh  of  the 
bias  6ub  changes  for  each  duty  cycle  k  according  to  the  following  definition: 


6ub  ~  j 

where  cr®  is  the  bias  uncertainty.  Since  only  the  variance  CTb  of  the  sensor’s  residual 
random  bias  is  specified,  the  bias  estimate  is  initialized  with  zero  mean  so  that  in  duty 
cycle  k  =  I,  6ub  ~  Note  that  these  error  models  would  apply  to  both  the 

acceleration  and  angular  rate  measurements;  for  simplicity,  this  discussion  is  limited  to 
the  accelerometer  measurements  along  a  single  axis.  The  bias  uncertainty  specifications 
for  the  conventional  INS  accelerometers  (cr^,^)  and  gyros  (cr^g)  are  calibrated  in  Appendix 
B  to  correspond  to  a  1  (km/hr)  conventional  navigation  quality  INS  position  error. 

Equations  (3.4)  and  (3.5)  are  combined  to  form  a  linear  regression  in  the  parameter  0, 
yielding  the  measurement  vector  z: 


z  =  H0  +  V 


(3.6) 


where  0  is  an  n  +  1  vector  of  true  acceleration  samples  and  the  accelerometer’s  residual 
bias  error  5ub  as  follows: 

^N-n 


0 


A 


6ub 


(n+2)xl 


25 


The  associated  regressor  matrix,  H,  is  given  as 


Ifi+l 

e 

H  = 

h" 

0 

Olx(n+l) 

1 

(n+3)x(«+2) 


where  I„+i  is  the  (n+l)x(n+l)  identity  matrix,  eisa(n+l)xl  vector  of  ones,  and  the 
transpose  of  the  vector  h  is 


h"  = 


1 

2 


1  ...  1 


1 

2 


lx(n+l) 


The  measurement  vector  in  the  linear  regression,  z,  consists  of  the  n  +  1  conventional 
INS  measurements,  the  single  CAI  accelerometer  measurement,  and  the  bias  error  estimate 
available  at  the  beginning  of  the  current  duty  cycle,  e.g.  duty  cycle  k.  z  is  expressed  as 
follows: 


A 

z  = 


u 

nu 

6ub 


(INS) 


rriN 

(CAI) 

m 

(k-1) 


(«+3)xl 


The  equation  error  vector  in  the  linear  regression,  v,  is  given  by 


6u 


(INS) 

N-n 


6u 


(INS) 

N 


Jk-i) 

_ 


-*(n+3)xl 


(3.7) 
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where  the  residual  variable  ~  A/'^O,  j.  Thus  equation  (3.6)  has  a  diagonal 

covariance  matrix  R,  expressed  as  follows  [14]: 


R  =  E  [vv^] 


cP-  T 

1  0(„+l)xl 

-1 - 

0(„+i)xl 

0lx(n+l) 

2  2 

«  (^CA/ 

0 

0lx(n+l) 

0 

(3.8) 


■'(«+3)x(o+3) 

Solving  the  linear  regression  in  (3.6)  according  to  the  following  equation  for  the  optimal 
estimate  vector  0: 

0  =  (H  R  H)  ‘  H^R  z  (3.9) 


yields  the  optimal  estimates  un-h,  ...,  mw  of  the  n  + 1  acceleration  samples  during  the  A:th  CAI 

- — -(.k) 

duty  cycle,  u{{k  -  a)AT  +  i6t),  i  =  0, ...,  n  and  the  accelerometer’s  bias  error  estimate  5ub 
obtained  at  the  end  of  the  kth.  CAI  duty  cycle.  The  covariance  of  the  parameter  estimation 
error  is  [14]: 


P  =  E  [(0  -  0)(0  -  0)^] 

=  (ffR-‘H)  ‘  (3.10) 


recalling  that  0  from  (3.9)  is: 

—  A-  —  - — ■(k)V 

&  =  \Ubl-n,-,UN,5Ub  I  (3.11) 

Since  the  inputs  u  were  treated  as  scalars,  0  will  need  to  computed  six  times  for  each  of 
the  three  acceleration  and  gyro  measurements.  However,  the  covariance  matrix  P  from  the 
linear  regression  is  the  same  for  all  acceleration  components.  Likewise,  P  will  be  the  same 
for  all  three  angular  rate  components.  Note  that  the  bias  uncertainty  featured  in 

equation  (3.8)  is  determined  from  P  as  follows: 

[erf)  =  P„+2,n+2,  k  =  2, ...;  =  cr^ 
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and  is  initialized  at  A:  =  1  from  the  respective  accelerometer  and  gyroscope  bias 
specifications,  and  ct^q,  both  derived  in  Appendix  A.  To  calculate  P,  substitute  for 
the  H  and  R  matrices  in  (3.10): 


H^R  H  = 


:  h  :  0, 


-n+l  I  I  V(„+i)xi 

Aoi'i' 


I  yj  I 


^_T 

_2  -■-n+l 

_/AS _ 

1 

0(„+l)xl 

0lx(n+l) 

1  1 

1  m2  ^2 

^CAl 

1 

0 

0lx(n+l) 

'  0 

1 

1 

Iii+i  I  c 

h"  :o 


0 


lx(/i+l)  I 


!  h  :  0, 

:  0  :  1 


(n+l)xl 


J_ 

INS_ 


.2  In+1 


-^e 

INS_ 


1 

"  ?ca;_  _ 


0 


0 

1 


lx(n+l) 


K-'T 


^_T  + _ 1 _ hh^  '  e 

^2  ^n+l  ^  ^2^2  ****  I  _2  ^ 

_ liVl 


9-2 

WAS 


n+l  _|_  1 


^2 

WAS 


K-‘T 


(n-i-2)x(n+2) 


(3.12) 


Since  (3.12)  is  in  block  form,  use 


Lemma  1:  Assume  A  is  symmetric  and  invertible.  The  inverse  of  the  blocked  matrix 

A  b 

b"  c 

where  the  scalar  d  =  c  -  b^A“^b.  □ 


A-^  +  U-^bb^A-* 

a 

1 

< 

-H|-93 

1 

-ib^A-i 

1 

d 

(3.13) 


28 


From  (3.12): 


A  =  hh^ 

o-Ins  o-f. 

b  =  e 

(t: 


INS 

n  +  I  1 


c  = 


+ 


cr 


n+l  1  1 

^  “  .2  ■*■  /  rt-n\2  “ 


cr;: 


■e'A-‘e 


0-7, 


To  solve  for  A  \  use  the  following  Matrix  Inversion  Lemma  (MIL): 


MIL:  Assume  the  dimensions  of  the  relevant  matriees  are  eompatible  and  the  required 
matrix  inverses  exist.  Then, 


(A.  -A,A;%y  =  A)‘  +A;^AM4-A,A;^A,yA,A;' 


□ 


Set 


A I  ^  I«+i 

^INS 

A2  =  h 
A3  =h^ 

A4  —  —n 


whieh  gives  the  inverse  of  A: 


A  cr^^^I„4ih(n  ^cai  ^  cr^^^I^+ih)  h 

erf 


/JVS  "+l  CAI 

0-4 

_2  T  .hh^ 


^;jvs^n+l 


^  cr,-.^;  +  erfv^h  h 
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Substituting  the  above  expressions  gives  the  covariance  of  the  parameter’s  estimation 
error  in  (3.10). 

Finally,  inserting  the  calculated  covariance  of  the  parameter’s  estimation  error,  P,  into 
equation  (3.9)  gives  the  accelerations’  estimate  and  current  accelerometer’s  bias  estimate: 


0 


I  h  ,  0,, 


■-n+1  I  I  '^(n+l)xl 


0 


Substitution  gives 


1  T 

2  An+1 

_INS _ 

1  ^(n+l)xl 

0(„+l)xl 

Olx(n+l) 

1  1 

1 

^CAI 

1 

0 

Olx(n+l) 

'  0 

1 

1 

0 


h  [a-‘  +  f  (A-‘bb^A-'  -  A-‘be")l 

rx/c  L  “  J 


^  INS 


h-  (b^A-'  -  e") 


^(A-'  +  lA-WA-)h 


-1  kta-1 


dna-l^, 


b^A  h 


dicr 


.(k-iy 


,A  b 


d(^c 


■n 


where 


A-Iu  u  ^INS 

b  =  cr^„,b  -  — ^ 

,  +  0-2  h  h 


hh"b 


CAI  '  INS' 
2 


ncr 


To  solve  for  first  expand  d 


n  +  \ 


+ 


=  e 


1 


+  0-^Ns(n  -  k) 


1 


cr 


INS 

1 


K"’) 


K-’) 

Then  invert  d  as  follows: 


+ 


cr 

n 


INS 

2 


(n  +  1)0-1,  - 


2  4 

^Ks 


n^cT^  (T^  h^h 

^  CAl  ^  ^  INS 


1 

d 


rdal^i  +  cr2„,  (n  -  i) 

+  («  -  y  ^Is) 


Concerning  the  recursion  for  cr\  : 


rdal^i  + 


ik). 


K‘‘‘f  +  {"  -  i) 


K’)’  =  ^ 


err, 


(3.14) 
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Then  from  equation  (3.14) 


(k) 

(tI’  = 


A 


n^o-^cA,  + 


(x-d 


err, 


n^crl^i  + 


cr 


ik-l) 


err 


Evidently,  cr®  is  monotonically  decreasing. 

Having  solved  the  linear  regression  in  (3.9)  to  obtain  the  estimates 
scalar  sensor  inputs  m„,  ...,  un  and  the  current  estimate  of  the  input  bias 


determined  by  substituting  (3.6)  into  (3.9): 


Un,  ...,  Mv  of  the 
,  the  error  6ui  is 


0  -  0  =  (H  R  H)  H  R  v 

Since  the  first  n  +  I  elements  of  0  6  consists  of  ...un,  it  is  now  known  that 

6ui  =  Hi  -  Ui 

=  (ff  R-‘H)-'  R-'v,  i  =  n, ...,  N 

where  ey  is  the  n  +  2  vector,  all  of  whose  elements  are  zero  except  the  y-th  element,  which 
is  1.  Recall  that  the  elements  of  the  n  +  3  equation  error  vector  v  in  (3.7)  consist  of  the 
errors  in  the  n  +  2  measurements  taken  during  the  CAI  measurement  part  of  the  duty  cycle 
and  the  last  element  is  the  current  knowledge  of  the  residual  uncertainty  (computed 
at  the  end  of  the  k  -  \  duty  cycle)  of  the  bias  error  6ub. 

Indeed,  at  the  beginning  of  the  current,  kth  duty  cycle,  information  concerning  the  bias 
error  is  supplied  as  follows: 

6ub  ~  N{6Ub  ) 

Furthermore,  the  statistics  of  are  expressed  as  follows: 

(3.15) 

Hence,  6ui  =  cr,  •  +  [a  linear  combination  of  the  first  u  +  2  elements  of  v]  where 

the  coefficients  or;  are  expressed  as  follows: 

cr,  =  e,,,_„(ffR-‘H)-‘ffR-'e„,3,  i  =  n,...,N  (3.16) 
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3.3  Proof  of  Concept 

The  combined  acceleration  and  angular  rate  estimates,  0,  during  the  CAI  measure¬ 
ment  interval  {{k-  \)lsT,klsT]  and  the  conventional  INS  acceleration  and  angular  rate 
measurements  prior  to  the  interval  are  used  in  an  inertial  frame  INS  to  calculate  the  ve¬ 
locity  and  position  x  of  the  aircraft.  This  is  accomplished  by  integrating  according  to  the 
simple  kinematic  differential  equation: 


X  =  u 

where  u  is  the  measured  acceleration.  The  angular  rotation  9  is  calculated  by  integrating 
the  angular  rate  to  as  follows: 

6  =  to 

The  error  in  position  is  caused  by  the  accelerometer  bias  and  any  platform 
misalignment  that  may  be  present.  The  misalignment  of  the  platform  due  to  the  gyroscopes’ 
errors  is  caused  by  the  gyro  bias  bg. 

The  position  and  velocity  estimation  errors  are  determined  by  the  conventional  and 
CAI  accelerometer  drifts  (cr,^5  and  trcAi,  respectively),  and  the  bias  b^  of  the  conventional 
accelerometer.  Specifically,  the  drift  induced  error  prior  to  the  CAI  accelerometer 
measurement  duty  cycle  is  determined  by  cr,„x,  and  the  drift  induced  error  during  the  CAI 
duty  cycle  interval  is  determined  from  the  covariance  matrix  P  in  equation  (3.10).  Similarly, 
the  angular  rotation  errors  are  determined  from  the  corresponding  INS  and  CAI  gyro  drifts 
and  gyro  bias  bo.  The  (n  +  2)x{n  +  2)  covariance  matrix  P  is  not  diagonal.  Further  analysis 
is  required  to  correctly  account  for  the  correlation  terms  in  P  and  solve  for  the  variance  of 
the  navigation  state  error  caused  by  the  acceleration  and  gyroscope  measurement  errors. 
For  both  the  INS  and  CAI  IMUs,  the  accelerometer  error  equations  are  as  follows: 

6x  =  6u,  6u  =  0 
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The  gyroscope  error  equation  is: 


66  =  6(j0 


3.4  Navigation 

A  wings  level  flight  at  constant  altitude  is  considered.  For  proof  of  concept,  the 
Earth  is  assumed  flat  and  non-rotating  and  the  inertial  navigation  frame  is  attached  to  the 
Earth.  The  simplified  discrete-time  navigation  state  error  equations  are  derived  from  their 
continuous-time  form  as  follows  [17]: 


6Xcit)  =  -t 

+  r. 

bG 

O-eMcit) 

where  the  navigation  state  is 


Xc 


|T 


x,y,z,  Vx,Vy,Vz,(p,  e,il/ 


(3.17) 


for  aircraft  position  coordinates  x,  y,  z;  velocities  v^,  Vy,  v^,  and  Euler  angles  6,  if/.  The 
continuous  time  state  space  matrices  and  are  given  as: 


0  I  0 

0  0 

II 

OOF 

II 

0 

0  0  0 

9x9 

_0  -c^ 

The  matrix  F  =  fx  is  the  skew  symmetric  matrix  form  of  the  nominal  specific  force  vector 
f.  During  wings  level  flight,  f  is  as  follows: 


Ux 


f  = 


0 

g 
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where  is  the  aircraft’s  longitudinal  acceleration  and  g  is  the  acceleration  of  gravity. 
Consequently,  the  matrix  F  is  as  follows: 

0  -g  0 

F  =  fx  =  ^  0  -u, 

0  Ux  0 

is  the  3x3  Direction  Cosine  Matrix  (DCM)  and  for  small  Euler  angles  is  as  follows: 

1  -(A  e 
^  1  -0 
-e  (p  1 

For  constant  altitude,  wings  level  flight  in  the  x  axis  direction,  C^  is  nominally  the  identity 
matrix,  I.  Pa  and  Pg  are  3x1  vectors  of  independent,  unity  Brownian  motions,  that  is, 
PA^t)  ~  Af  (0,  t  •  I)  and Pcit)  ~  Af  (0,  t  •  I).  The  rate  gyro  drift  parameter  cTg^  is  derived  from 
the  rate  gyro’s  drift  specification  in  Appendix  A. 

Converting  the  dynamics  in  (3.17)  to  discrete  time,  the  navigation  state  discrete-time 


error  equation  is 

kl 

=  Md6xi  +  Yd  -I-  Fc 

be 

where 

I  lAT  -fFAT^ 
Mrf  =  =0  I  -FAT 

0  0  I 


(3.18) 
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and  w^.  ~  N  (0,1),  Wc,  ~  A/"  (0,1).  The  discrete-time  accelerometer  and  gyro  parameter 
featured  in  (3.18),  and  respectively,  are  derived  in  Appendix  A. 

It  is  herein  assumed  that  intially  the  INS  is  perfectly  aligned  (no  alignment  error)  and 
therefore  the  navigation  state  error  is  exclusively  brought  about  by  the  accelerometers’  and 
gyros’  biases  and  drifts.  As  such,  the  navigation  state  error’s  statistics  will  be  determined 
by  the  accelerometers’  and  gyros’  bias  statistics  and  the  intensity  of  their  respective  drifts. 

After  the  measurement  part  of  the  kth  duty  cycle  is  complete  and  the  linear  regression 
in  (3.9)  has  been  solved,  one  reruns  the  mechanization  equations  which  are  modified  as 
follows: 

=  fd  (x^,,  u,) ,  i  =  n,  ...,N-  1  (3.19) 

where  u,  is  a  6  x  1  vector  of  the  scalar  accelerometer  and  gyro  estimates  Hn,  ...,'un-i 
obtained  from  (3.9).  The  result  is  a  smoothed  navigation  state  estimate  ...,Xm„  for 
the  CAI  measurement  portion  of  the  duty  cycle.  The  navigation  state  filtered  estimate 
is  obtained  in  real  time  with  some  computational  delay;  however,  the  smoothed  navigation 
state  estimates  x^^  j,  are  delayed  by  the  computational  delay  -i-  6t,  ...+,(n  -  l)6t, 

where  6t  are  the  discrete-time  instants.  In  view  of  the  mechanization  equation  expressed 
in  (3.19),  the  discrete-time  navigation  state  error  for  the  smoothed  navigation  state 
during  the  CAI  measurement  part  of  the  duty  cycle  is  expressed  as  follows: 

=  Mrfdx^,  -t  TdSvii,  i  =  n,  ...,N  -  1  (3.20) 
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where 


S^nij  —  ^nii  ^  ^  “1“  1 ,  . . . , 

6Ui  =  u,  -  u,-,  i  =  n, ...,  N  -  I 

and  dXm;  and  du,  are  6  x  1  vectors  of  the  scalar  accelerometer  and  gyro  measurement  errors. 
The  statistics  of  the  scalar  input  errors  6un,  ...,6un-i  are  provided  from  the  solution  to  the 
linear  regression  in  (3.10)  obtained  at  the  completion  of  the  kih.  duty  cycle,  and  are  as 
follows: 

Stifi 

~Af(0,P®) 

6un 
Ak) 

The  covariance  of  the  smoothed  navigation  state  estimation  error  expressed  in 
(3.20),  may  be  rewritten  as  follows: 

j-n 

=  Mf  ^  Mf  ”“'rddu„+/_i,  j  =  n  +  l, ..., N  (3.21) 

i=i 

where  dXm„  ~  A/'(0,P|f'^^j.  The  covariance  of  the  smoothed  navigation  state  estimation 
error  ...,  P^^^  is  then  expressed  as: 

P“'>  S  E  +  r,p®rl  +  M,E  («x„>u:)  Fj 

+  rj(E((Sx„/uI)f  Mj 

Pf’  S  E(ix„<,) 

To  solve  for  the  covariances  when  i  =  n,...,N  -  1,  dx„,^  must  first  be 

considered.  For  the  first  half  of  each  duty  cycle  (prior  to  the  CAI  measurement  portion). 
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the  mechanization  equation  for  the  CAI-aided  INS  is  expressed  as  follows: 

Am,+i  —  Jdy^mn^mi  ) 

The  discrete  time  navigation  state  error  equation  is: 


S'^nii+i  -  +  T +  r, 


O’v^Wa,- 


,  i  =  0,  ...,n  -  1 


where  ^Ms  a  6  x  1  vector  of  the  scalar  residual  uncertainties,  defined  by  (3.15),  of 
the  accelerometer  and  gyro  bias  errors.  Note  that  Thus 


6x,n^  =  + 


'  n-l  ' 

V  1=0  ^ 

The  correlations  E  ^dx„du^j  for  i  =  n, ...,  N  -  I  can  then  be  calculated  as  follows: 


/n-l  \ 


E(dx„du^)  =  gM'Jr,E(^-'^a4^-i^) 

/  n-\ 


V  1=0  ) 


where  cr,-  was  given  in  (3.16). 

Returning  to  (3.18),  if  it  is  assumed  that  the  accelerometer  and  gyro  errors  and  be 
respectively  are  constant,  random  Gaussian  distributed  biases,  the  state  error  vector  dx  may 
be  augmented  with  the  vectors  b^  and  b^  as  follows: 


dx; 


6xi,  6yi,  6zi,  Svjc„  6vy.,  6v^.,  6(pi,  69i,  di/f,-,  b^,.,  be, 


Consequently,  the  dynamics  matrix  is  augmented  with  discrete  time  matrix  as  follows: 


A  = 


T, 

^6x9  lexe 
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and  Fc  is  augmented  to  acknowledge  the  additional  bias  states  to  form  discrete  time  matrix 
G: 


06x6 


0  0 
I  0 
0  I 
0  0 
0  0 


15x6 


The  final  discrete-time  form  of  the  augmented  navigation  state  error  dynamics  is  expressed 
as  follows: 

dx,+i  =  A6xi  +  G  VQw,- 


where 


Vq 


0 

Wa,- 

,W;  = 

0 

1— 1 

b 

6x6 

Wg, 

3.5  Simulation 

For  simulation  of  the  conventional  INS,  (3.18)  is  modified  slightly  to  reflect  the 
calibration  performed  in  Appendix  B  as  follows: 


dXi+i  =  MddXi  -t  Frf 

b. 

+  r^ 

bg 

(3.22) 


where  the  accelerometer  and  gyro  errors  and  bg  are  3  x  1  vectors  of  constant,  random 
Gaussian  distributed  bias,  so  b^  ~  Af^O,  and  bg  ~  Af^O,  cr^^lj.  The  values  for  the 
bias  variances,  and  cr^^,  and  discrete-time  accelerometer  and  gyro  parameters,  and 
cTg^  respectively,  are  calculated  in  Appendix  B. 
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The  CAI-aided  INS  simulation  mechanization  equations  for  the  first  N  -  n  measure¬ 
ments  of  each  duty  cycle  k  =  1, 3600  are  as  follows: 

-t  Td  -  6\xl  j ,  5x„„  =  0  6  i  =  0, N-n-l  (3.23) 
where  is  a  6  x  1  vector  consisting  of  each  scalar  conventional  INS  accelerometer  and 

-  (k—i) 

gyro  measurement  described  in  (3.5).  Likewise,  is  a  6  x  1  vector  of  each  of  the 

scalar  accelerometer  and  gyro  bias  error  estimates.  Note  that  when  k  =  1,  all  six  bias  error 
estimates  are  initialized  as  zero. 

The  CAI-aided  INS  simulation  mechanization  equation  for  the  CAI-measurement 
portion  of  each  duty  cycle  k  is  as  follows: 

dx^,^,  =  +  TdUi,  i  =  N-n, ...,  N  (3.24) 

where  u,  is  a  6  x  1  vector  of  the  estimated  scalar  accelerometer  and  gyro  measurements 
provided  from  0  in  the  linear  regression  expressed  in  (3.9). 

3.6  Revision  of  Model  1 

Vanderbruggen  and  Mitchell  [18]  have  recently  proposed  a  numerical  method  to 
generate  a  collimated,  continuous  source  of  pin-polarized  atoms  which  could  be  used  for 
cold  atom  interferometry.  Increasing  the  rate  of  available  atoms  should  likewise  increase 
the  measurement  availability.  With  this  concept  in  mind,  the  initial  integration  model  is 
revised  to  investigate  what  improvements,  if  any,  would  be  gained  from  increasing  the  CAI 
measurement  rate  while  maintaining  the  50  percent  duty  cycle.  CAI  measurements  were 
simulated  twice  per  second  and  reflected  the  average  of  the  acceleration  and  angular  rate 
measurements  during  the  two  measurement  duty  cycles,  as  illustrated  in  Figure  3.2.  Finally, 
the  model  is  revised  once  more  to  simulate  a  5  Hz  CAI  measurement  rate  with  50  percent 
duty  cycle.  All  model  results  are  collected  and  tabulated  in  Chapter  4. 
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Figure  3.2:  Revised  Model  1.  Acceleration  over  two  CAI  duty  cycles. 


3.7  Alternate  Model 

An  alternate  model  is  considered  in  an  attempt  to  improve  the  CAI-aiding  performance 
inspired  from  research  conducted  by  Omr  et  al  [9] .  As  discussed  in  Section  I,  Omr,  Georgy 
and  Noreldin  integrate  position  and  velocity  measurements  provided  from  multiple  MEMs 
INS.  The  MEMs  devices  calculate  ownship  position  and  velocity  by  integrating  their  unique 
acceleration  measurements.  The  alternate  model  assumes  that  the  CAI  system  operates 
as  an  inertial  measurement  unit  with  its  processor  for  calculating  position,  velocity  and 
angular  rate  states,  Eikewise,  the  conventional  INS  unit  has  its  own  integration 

processor  which  outputs  position,  velocity  and  angular  rate  states  The  conventional 

INS  and  CAI  IMU  state  measurements  are  optimally  combined  using  basic  Kalman  gain 
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principles  (discussed  in  Chapter  2)  to  calculate  the  best  estimate  position,  velocity,  and 
angular  rate  states  x.  The  best  estimate  x  is  then  substituted  for  the  previously  calculated 
x(iNS)  in  the  conventional  INS  mechanization  equations,  as  shown  in  Figure  3.3. 


Figure  3.3:  Model  2  System  Concept. 


The  9x1  INS  and  CAI  state  measurement  vectors  are  defined  as  follows: 

x(INS)  ^  ^ 

x(CAI)  _  x(TRUE)  y 
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where  is  the  true  measurement  state  vector,  and  the  INS  and  CAI  measurement 

error  vectors  are  ^  and  v,  respectively.  The  errors  are  assumed  to  be  zero-mean,  stochastic 
variables  with  covariances  P,  as  expressed  with  the  following  equations: 

V  ~  yv(0, 

The  INS  and  CAI  state  measurement  vectors  may  be  augmented  as  z  as  follows: 

z  =  Hx  -I-  v 

x(INS) 

x(CAI) 

*-  -*18x1 

with  regressor  matrix  H  and  augmented  error  vector  v 

I9x9  ^ 

H=  ,  v  = 

I9x9  J  [ 

and  augmented  covariance  matrix  R 

p(INS)  I  Q 
r  V9x9 

R=  t 

0  !  p(CAI) 

V9x9  I  -T 

*-  I  -*18x1 

From  the  Kalman  filter  optimization  discussed  in  Chapter  2,  the  estimated  state  vector 

is  calculated  from  the  following: 

r  -1-1 

X  =  H  iRH  Hz 

where  K  is  the  gain,  defined  as: 

p(INS)jjT 

^  jjp(INS)jj  p(CAI) 

This  concept  model  would  work  if  the  CAI  measurements  were  consistently  available. 
However,  if  there  were  CAI  accelerometer  or  gyro  measurement  dropouts  due  to  high 
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dynamics,  the  model  is  no  longer  useful  beeause  the  CAI  meehanization  equations  would 
have  to  be  re-initialized.  For  this  reason,  this  alternate  model  was  not  pursued  beyond  a 
eoneeptual  design. 
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IV.  Results  and  Analysis 


4.1  Monte  Carlo  Analysis 

A  Monte  Carlo  (MC)  simulation  of  1000  runs  is  performed  to  validate  the  navigation 
algorithm’s  performanee  over  one  hour.  The  aireraft’s  position  and  veloeity  are  initialized 
at  0  meters  and  0  (m/s),  respeetively.  The  ealibrated  variances  of  the  measurement 
uncertainties  are  given  in  Table  4.1. 


Table  4.1:  Sensor  Uncertainties  for  Scenario  1  =  1000  (m/hr),  Scenario  2  =  1500  (m/hr) 


Scenario  I 

Bias  Variance 

Measurement  Noise  Variance 

INS  accelerometer 

(7.72  X  10-5  m/s^)^ 

(l.27x  10-2  m/s2/Viu)^ 

INS  gyro 

(6.55  X  10"^  rad/hrj 

(9.27  X  10-2  rad/Vhr)^ 

CAI  accelerometer 

N/A 

(2.00  X  10-5  m/s2/  Vni)^ 

CAI  gyro 

N/A 

(l.47  X  I0-‘’  rad/Vhr)^ 

Scenario  2 

INS  accelerometer 

(l.I6x  10-^  m/s^)^ 

(l.90x  10-2  m/s2/Vlu)^ 

INS  gyro 

(9.83  X  I0-‘'  rad/hr)^ 

(l.39x  lO-'^  rad/Vhr)^ 

CAI  accelerometer 

N/A 

(2.00  X  10-5  m/s2/  Viu)^ 

CAI  gyro 

N/A 

(l.47  X  I0-‘’  rad/Vhr)^ 

From  this  analysis,  it  is  determined  that  aiding  a  conventional  INS  with  CAI-based 
acceleration  measurements,  or  vice  versa,  aiding  a  CAI-based  INS  with  conventional  INS 
measurements,  reduced  the  navigation  error  by  approximately  48.5  percent,  as  illustrated 
in  Figure  4.1.  An  additional  Monte  Carlo  simulation  is  run  with  the  conventional  INS 
position  accuracy  calibrated  to  1500  (m/hr),  reflecting  navigation-grade  INS  specifications. 
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to  see  if  the  CAI  aiding  would  be  as  effeetive.  As  shown  in  Figure  4.2,  the  standard 
deviation  is  reduced  by  approximately  47.2  percent,  which  suggests  that  the  unaided  system 
performance  does  not  significantly  impact  the  CAI-aiding  capabilities  when  used  with 
the  proposed  mechanization  algorithm.  Table  4.2  shows  the  unaided,  aided,  and  percent 
correction  for  the  position,  velocity,  and  attitude  parameters  in  both  scenarios. 


Figure  4.1:  Monte  Carlo  Simulation  results  for  Position  Error  over  1000  runs.  Calibrated 
standard  deviation  position  error  for  conventional  INS  is  1000  (m/hr). 


4.1.1  Varying  Conventional  INS  Data  Rate. 

The  original  simulation  assumes  that  the  CAI  measurements  would  be  aiding  a 
conventional  INS  with  10  Hz  data  rate.  This  assumption  is  made  to  simplify  the  algorithm 
development;  however,  conventional  INS  systems  operate  at  much  higher  data  rates  in 
the  field,  typically  in  the  range  of  50  to  300  Hz.  To  account  for  this  disparity,  the  CAI- 
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Figure  4.2:  Monte  Carlo  Simulation  results  for  Position  Error  over  1000  runs.  Calibrated 
standard  deviation  position  error  for  eonventional  INS  was  1500  (m/hr) 


aiding  algorithm  is  tested  with  eonventional  INS  sampling  rates  of  50  Hz,  100  Hz  and  150 
Hz.  Computer  memory  limitations  prevented  simulations  above  150  Hz.  The  INS  drift 
and  bias  uneertainties  are  re-ealibrated  for  1500  (m/hr)  drift  to  eompensate  for  eaeh  of 
the  differing  data  rates  as  shown  in  Table  4.3.  Note  that  the  CAI  measurement  data  rate 
is  kept  at  1  Hz,  so  no  reealibration  was  neeessary.  The  unaided  and  aided  RMS  values 
for  the  position,  veloeity  and  attitude  are  given  in  Table  4.4.  The  original  10  Hz  results 
are  also  ineluded  for  eomparison  purposes.  Figure  4.3  illustrates  an  inverse  relationship 
between  the  error  eorreetion  pereentages  and  eonventional  INS  data  rates.  As  the  data  rates 
inerease,  the  INS  measurements  are  weighted  more  heavily  and  overshadow  the  single  CAI 
measurement.  With  more  samples,  there  is  more  potential  for  anomolous  noise  to  enter  the 
INS  measurements  whieh  may  not  be  aeeounted  for  by  the  algorithm’s  weighting  matrix  R. 
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Table  4.2:  Unaided,  Aided  Error  and  Error  Correction  Values  after  One  Hour  for  Varied 
INS  Drifts 


1000  (m/hr)  INS  drift 

Position  (m) 

Velocity  (m/s) 

Attitude  (rad) 

Unaided  INS  (RMS) 

1051.25 

0.68 

2.98  X  10-5 

CAI-Aided  INS  (RMS) 

541.85 

0.34 

1.40  X  10-5 

Error  Correction 

48.46  % 

50.19% 

53.19  % 

1500  (m/hr)  INS  drift 

Unaided  INS  (RMS) 

1554.50 

1.00 

4.60  X  10-5 

CAI-Aided  INS  (RMS) 

819.60 

0.50 

2.22  X  10-5 

Error  Correction 

47.23  % 

50.18  % 

51.70% 

However,  the  drop  in  aiding  performance  is  marginal;  with  a  fifteen-fold  increase  in  data 
rate,  the  error  correction  percentage  decreases  by  approximately  six  percent. 


Table  4.3:  Calibrated  Conventional  INS  Uncertainties  for  Different  Data  Rates 


50  Hz 

Bias  Variance 

Measurement  Noise  Variance 

INS  accelerometer 

(l.I6x  10-^ 

(4.25  X  10-2  m/s2/ 

INS  gyro 

(9.83  X  I0-‘^  rad/hr)^ 

(3.1 1  X  10-*’  rad/Vhr)^ 

100  Hz 

INS  accelerometer 

(l.I6  X  10  m/s^j 

(6.01  X  10-2  m/s2/  Viu)^ 

INS  gyro 

(9.83  X  I0-‘^  rad/hr)^ 

(4.39  X  10-*’  rad/Vhr)^ 

150  Hz 

INS  accelerometer 

(l.I6x  10-^  m/s^)^ 

2 

(7.37  X  10-2  m/s2/  Viu) 

INS  gyro 

(9.83  X  I0-‘^  rad/hr)^ 

(5.39  X  10-*’  rad/Vhr)^ 
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Figure  4.3:  Position  Error  Corrections  from  Monte  Carlo  Simulations  for  Varying 
Conventional  INS  Data  Rates. 


4.1.2  Varying  Acceleration  Input. 

In  the  field,  the  CAI-aided  INS  will  be  subjected  to  varying  accelerations.  To  test  its 
performance,  two  scenarios  are  simulated  where  the  accelerations  are  0.5  Hz  sinusoidal 
inputs.  The  first  acceleration  waves  are  given  an  amplitude  of  1.5  g’s  and  the  second,  3 
g’s.  The  3  g  amplitude  represents  the  maximum  cut-off  for  low  dynamics,  based  upon  the 
author’s  field  test  experience  and  previous  research  from  Canciani  [12].  CAI  measurements 
are  flagged  as  unusable  when  subjected  to  accelerations  above  3  g’s  because  the  vibrations 
saturate  the  sensitive  instruments  of  the  low-bandwidth  system.  It  should  be  noted  that 
flagging  an  aiding  system’s  measurements  in  poor  environmental  conditions,  e.g.  poor 
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Table  4.4:  Unaided  and  Aided  RMS  Error  after  One  Hour  for  Varied  INS  Data  Rates 


10  Hz 

Position  (m) 

Veloeity  (m/s) 

Attitude  (rad) 

Unaided  INS  (RMS) 

1554.50 

1.00 

4.60  X  10-5 

CAI-Aided  INS  (RMS) 

819.60 

0.50 

2.22  X  10-5 

Error  Correetion 

47.23  % 

50.18  % 

51.70% 

50  Hz 

Unaided  INS  (RMS) 

1529.13 

1.02 

4.30  X  10-5 

CAI-Aided  INS  (RMS) 

876.21 

0.56 

2.21  X  10-5 

Error  Correetion 

42.70  % 

45.61  % 

48.55  % 

100  Hz 

Unaided  INS  (RMS) 

1566.95 

1.03 

4.52  X  10-5 

CAI-Aided  INS  (RMS) 

897.16 

0.57 

2.34  X  10-5 

Error  Correetion 

42.74  % 

44.53  % 

48.20  % 

150  Hz 

Unaided  INS 

1505.09 

0.99 

4.46  X  10-5 

CAI- Aided  INS 

889.55 

0.56 

2.35  X  10-5 

Error  Correetion 

40.90  % 

43.38  % 

47.42  % 

GPS  GDOP,  is  a  eommon  praetiee  in  industry.  For  eaeh  seenario,  the  INS  and  CAI  drift 
and  bias  uneertainties  are  re-ealibrated  for  1500  (m/hr)  and  5  (m/hr)  respeetively,  beeause 
the  values  in  the  transition  matrix  ehanged  for  different  aeeeleration  inputs.  The  CAI  data 
rate  are  limited  to  1  Hz  while  the  INS  data  rates  are  limited  to  10  Hz.  The  unaided,  aided, 
and  error  eorreetion  values  for  the  position,  veloeity  and  attitude  are  given  in  Table  4.5. 
Comparing  the  two  seenarios,  the  CAI-aided  INS  has  slightly  more  aeeurate  performanee 
when  the  aeeeleration  amplitude  is  inereased  to  3  g’s.  This  is  possibly  due  to  the  CAI-aiding 
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being  more  effective  when  larger  conventional  INS  errors  are  introduced  from  increased 
dynamics.  The  growth  of  INS  errors  is  highly  dependent  upon  the  vehicle  trajectory,  as 
seen  in  the  time-dependent  transition  matrix  A.  More  scenarios  should  be  tested  before 
making  conclusions  regarding  the  acceleration  impacts.  It  should  also  be  emphasized  that 
this  trend  would  likely  not  hold  in  accelerations  greater  than  3  g’s.  Multiple  error  sources 
are  not  modeled  in  this  research  that  would  become  more  influential  in  high  dynamics,  e.g. 
accelerometer  and  gyroscope  scale  factors.  Additionally,  current  CAI  technology  limits  its 
measurement  bandwidth  to  dynamics  to  less  than  3  g’s. 


Table  4.5:  Unaided,  Aided  Error  and  Error  Correction  Values  after  One  Hour  for  Varied 
Acceleration 


1.5  g  Acceleration 

Position  (m) 

Velocity  (m/s) 

Attitude  (rad) 

Unaided  INS  (RMS) 

1579.15 

1.04 

4.76  X  10-5 

CAI-Aided  INS  (RMS) 

832.08 

0.52 

2.22  X  10-5 

Error  Correction 

47.31  % 

49.62  % 

53.53  % 

3  g  Acceleration 

Unaided  INS  (RMS) 

1623.58 

1.07 

4.62  X  10-5 

CAI-Aided  INS  (RMS) 

826.07 

0.52 

2.16  X  10-5 

Error  Correction 

49.12  % 

51.06% 

53.15  % 

4.1.3  Varying  CAI  Measurement  Data  Rate. 

As  discussed  in  Chapter  3,  recent  research  advances  in  cold  atom  physics  have  shown 
the  potential  for  having  continuously  sourced  atoms,  which  would  allow  for  higher  sample 
rates  and  bandwidth  [18].  To  test  the  impact  of  increased  CAI  sample  rates,  Monte 
Carlo  simulations  are  run  with  2  Hz  and  5  Hz  CAI  data  rates.  The  conventional  INS 
sample  rates  are  maintained  at  10  Hz.  The  unaided,  aided,  and  error  correction  values 
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for  the  position,  velocity  and  attitude  are  given  in  Table  4.6,  with  the  1  Hz  CAI  data 
rate  results  included  for  comparison  purposes.  A  five-fold  increase  in  CAI  data  rates 
gave  remarkable  improvements  in  navigation  accuracy,  including  a  nearly  87  percent 
correction  of  RMS  position  drift  error.  These  results  demonstrate  how  improvements  in 
CAI  measurement  availability  would  significantly  increase  the  effectiveness  of  CAI-based 
aiding  when  implementing  the  proposed  mechanization  algorithm. 


Table  4.6:  Unaided,  Aided  Error  and  Error  Correction  Values  after  One  Hour  for  Varied 
CAI  Data  Rates 


1  Hz  CAI 

Position  (m) 

Velocity  (m/s) 

Attitude  (rad) 

Unaided  INS  (RMS) 

1554.50 

1.00 

4.60  X  10-5 

CAI-Aided  INS  (RMS) 

819.60 

0.50 

2.22  X  10-5 

Error  Correction 

47.23  % 

50.18  % 

51.70% 

2  Hz  CAI 

Unaided  INS  (RMS) 

1577.96 

1.04 

4.75  X  10-5 

CAI-Aided  INS  (RMS) 

836.52 

0.53 

2.24  X  10-5 

Error  Correction 

49.99  % 

49.05  % 

52.91  % 

5  Hz  CAI 

Unaided  INS  (RMS) 

1533.84 

1.00 

4.34  X  10-5 

CAI-Aided  INS  (RMS) 

206.57 

0.22 

1.35  X  10-5 

Error  Correction 

86.53  % 

77.76  % 

68.96  % 
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V.  Conclusions  and  Recommendations 


5.1  Conclusion 

High  accuracy  cold  atom  interferometers  offer  a  potential  quantum  leap  in  inertial 
navigation  accuracy  and  possibly  an  autonomous  navigation  alternative  for  GPS  navigation. 
However,  the  CAI  accelerometers  and  gyros’  low  duty  cycle  significantly  reduces  its 
viability  as  a  standalone  INS.  Fusing  the  CAI-based  accelerometer  and  gyroscope 
measurements  together  with  conventional,  navigation  grade  INS  measurements  increases 
the  overall  system’s  bandwidth  while  also  reducing  the  navigation  error.  It  should  be  noted 
that  the  covariance  calculations  add  processing  delay,  with  a  minor  reduction  in  duty  cycle 
(dependent  upon  the  processor  capabilities).  In  this  research,  a  10  Hz  free-INS’s  position 
errors  are  reduced  by  47.2  percent  by  employing  CAI  measurements,  despite  the  one  second 
gaps  in  their  availability.  When  the  conventional  INS  data  rate  is  increased  to  150  Hz, 
the  CAI-aided  INS  showed  a  marginal  decrease  in  navigation  accuracy.  Increasing  the 
CAI  measurement  sampling  rate  made  a  significant  improvement  in  the  CAI-aided  INS 
performance. 

5.2  Recommendations  for  Future  Research 

Further  investigations  should  be  made  to  determine  the  CAI-aiding  performance  dur¬ 
ing  dynamic  flight  scenarios  where  the  CAI  may  incur  longer  gaps  between  measurements, 
similar  to  research  conducted  by  Canciani  and  Raquet  [11]. 

As  explained  Section  1.4,  the  acceleration  and  gyro  measurement  models  are 
simplified  to  include  a  single  bias  and  drift  induced  error.  The  Earth  is  also  assumed  to 
be  flat  and  non-rotating  for  the  navigation  simulations.  More  complex  error  and  mechanics 
models  are  proposed  in  [11]  and  [19].  These  models  could  be  employed  to  make  a  realistic 
evaluation  of  their  contributive  effects  in  more  dynamic  flight  scenarios. 


52 


Errors  caused  from  the  processing  delay  or  lag  should  be  investigated  in  a  real-time 
model.  The  errors  caused  from  gaps  in  acceleration  and  angular  rate  measurements  become 
more  significant  when  those  corresponding  dynamics  are  rapidly  changing,  particularly 
in  medium  to  high  dynamics.  The  developed  algorithm  performance  could  be  compared 
with  an  extended  Kalman  filter  which  estimate  the  INS  mechanization  during  gaps  in  CAI- 
aiding.  Increased  CAI  measurement  availability  may  reduce  the  impact  from  these  delays. 

Cold  atom  technology  is  developing  at  a  fast  pace  due  to  its  application  in  multiple 
areas  where  precision  measurements  are  a  necessity.  As  more  CAI  improvements  and 
demonstrations  are  made,  system  limitations  and  parameters  should  be  updated  to  reflect 
actual  performance.  Ideally,  flight  testing  of  actual  CAI  hardware  could  be  performed  in 
the  future  to  test  advanced  mechanization  models  and  their  aiding  algorithms. 
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Appendix  A:  Analytic  Sensor  Drift  Modeling 


A.l  Gyro  Drift 

The  continuous-time  rate  gyro  dynamics  are  modeled  as  a  stochastic  differential 
equation 

m  =  coit)  +  b^  +  (Tj{t),  0(0)  =  0  (A.l) 

where  m(0  is  the  true  input  angular  rate,  bo  is  the  gyro’s  residual  random  bias  bg  ~ 
N  (O,  (tI^^  and  the  stochastic  process  /3(-)  is  a  unit  Brownian  motion,  that  is,  /3{t)  ~  N  (0,  t). 
Note  that  the  units  of  the  strength  of  the  gyro’s  drift  parameter  cTg^  are  (rad/  Vs).  If  a)(t)  =  0 
and  bo  =  0,  the  solution  to  the  stochastic  differential  equation  (A.l)  is  as  follows: 

(A.l) 


The  specified  rate  gyro’s  measurement  uncertainty  caused  by  drift  cTg  is  defined  as  the 
angular  position  uncertainty  brought  about  by  integrating  the  rate  gyro’s  output  over 
one  hour.  In  industry,  it  is  commonly  referred  to  as  the  Angle  Random  Walk  (ARW) 
specification  because  the  gyro  drift  introduces  a  zero-mean  random  walk  error  in  the  angle 
[10].  In  other  words,  if  the  true  angular  rate  a)(t)  =  0,  then  at  one  hour,  the  following  is 
true: 

0(3600)  ~A^(0, 0-2)  (A.3) 

Relating  equations  (A.2)  and  (A.3)  when  t  =  3600,  and  using  the  manufacturer’s  provided 
gyro  specification  cTg,  we  obtain 


O-g 

60. 


(rad/  Vs) 


(A.4) 


which  defines  the  parameter  cTg^  featured  in  the  continuous-time  gyro  dynamics  modeled 
by  the  stochastic  differential  equation  in  (A.l)  and  in  the  continuous-time  INS  navigation 
state  error  equation  in  (3.17). 
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The  discrete-time  form  of  (A.l)  is  as  follows: 


=  Oi  +  (a>i  +  bc)AT  +  o-g^rii,  Oq  =  0,  i  =  0, 1, ... 


(A.5) 


where  the  intensity  of  the  process  noise,  cr^^,  is  the  discrete-time  equivalent  of  the  parameter 
(Tg^.  The  random  variable  ?/,  is  stipulated  as  ?7,  ~  AfCO,  1).  When  m,-  =  0  and  be  =  0,  the 
solution  to  the  discrete-time  stochastic  difference  equation  (A.5)  is: 


=  ^0,  Yj 

i=0 

~Af(0,zcr^J,  z  =  0, 1,... 


(A.6) 


Since  AT  is  related  to  the  sampling  rate  fs  through  the  following  equation: 


AT  = 


then  after  one  hour,  i  =  3600  •  fs,  where  we  recall  the  units  of  Hz  are  (1/s).  From  (A. 3), 
(A. 4),  and  (A.6),  the  parameter  featured  in  the  stochastic  difference  equation  (A.5)  is  as 


follows: 


•  Vat 


Assuming  the  gyro  drift  and  gyro  bias  uncertainties  are  the  only  error  sources,  the 
continuous-time  navigation  state  error  equation  (3.17)  would  be  specified  as  follows: 


6xc(t)  =  MgSXcit)  -H  -h  o-gTcMt) 


(A.7) 


where  ha  ~  N  (O,  cr^gl),  ;Sg  is  a  3  x  1  vector  of  independent  unity  Brownian  motions,  and 


Tg.  =  0 
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Note  that  in  (A. 7)  the  units  of  are  (rad/s).  If  is  specified  in  terms  of  (deg/hr),  then 
convert  to  (rad/s)  using  the  following  equation: 

(rad/s)  :=  ^ 

7T 

The  discrete-time  form  of  the  navigation  state  error  equation  (A.7)  is 

dXi+i  =  M^dXi  +  re^be_  -t  O-.J'a^y/cr  ^*0  =  0  (A.8) 

where 

Wg,  is  a  3  X  1  vector  of  independent,  Gaussian  random  variables  with  Wg,  ~  N  (0, 1)  for 
i  =  0, 1, ...  and  bg,  is  a  3  x  1  vector  of  constant,  random  biases,  that  is,  bo,^;  =  bo^,  bgo  ~ 
Af^O,  cr^^lj.  The  constant  gyro  biases  may  be  augmented  into  the  state  vector,  such  that 
(A.8)  becomes 


where 


dx,+i  =  Aedx;  -r  0-0^  GgW, 


dxn  ~  N 


0, 


0  0 


0  <Ijy 


Tg. 

Ag  = 

,  Go  -  CTg^ 

0 

I 

0 

A.2  Accelerometer  Drift 

The  continuous-time  accelerometer  dynamics  are  modeled  as 


(A.9) 


v(t)  =  u{t)  +  bj,  +  v(0)  =  0  (A.  10) 

where  u(t)  is  the  true  acceleration,  is  the  accelerometer’s  random  bias  ~  Af^O, 
and  yS(t)  is  a  unit  Brownian  motion.  Hence,  the  units  of  cr„^  are  (m/s  Vs).  If  the  true 
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acceleration  u{t)  =  0  and  bj^  =  0,  the  solution  of  the  stoehastie  differential  equation  in 
(A.  10)  is  as  follows: 

v(t)-N{0,cTlt)  (A.ll) 

The  aeeeleration  measurement  uneertainty  eaused  by  drift  cr„  (m/s)  is  speeified  as  the 
veloeity  uneertainty  brought  about  by  integrating  the  aeeelerometer’s  output  over  one  hour. 
The  aeeelerometer  drift  introduees  a  zero-mean  random  walk  error  in  the  veloeity  output, 
so  it  is  aptly  ealled  the  Veloeity  Random  Walk  (VRW)  speeifleation  [10].  In  other  words, 
if  the  true  aeeeleration  u{t)  =  0,  then  at  one  hour: 

v(3600)  ~  Af(0,(r2)  ^2) 

From  (A.ll)  and  (A.  12),  when  t  =  3600,  cr„^  is  as  follows: 

=  (m/sV5) 

thus  ealeulating  the  parameter  whieh  features  in  the  stoehastie  differential  equation 
(A.  10)  using  the  manufaeturer  supplied  speeifleation  cr„. 

Converting  the  dynamies  (A.  10)  to  diserete-time  form  gives  the  following  equation: 


Vi+i  =  Vi  -I-  (Ui  +  bA)AT  +  vq  =  0  (A.  13) 

where  is  the  diserete-time  equivalent  of  and  as  before,  ?/,  ~  A/'(0, 1).  When  m,  =  0 
and  b/,  =  0,  the  following  is  true: 

~  Af(0,/cry,  /  =  0,1,...  (A.14) 


From  (A.ll),  (A.  12),  and  (A.14),  the  parameter  featured  in  the  stoehastie  differenee 
equation  (A.  13)  is: 


cr 


(m/s) 


Finally,  ineluding  the  aeeelerometer  bias  and  drift  error  sourees  in  the  navigation  state 
error  dynamies  equation  (A.7)  gives  the  following  equation: 


dXc(0  -  McdXc(0  -I-  Fe^bG  +  cr -I-  cr„^rA^^A(0 


(A.15) 
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where 


0 

r.,  =  I 

0 

;Sa  is  a  3  X 1  vector  of  independent  unity  Brownian  motions,  and  b,i  ~  Af  (O,  Note  that 
the  accelerometer  bias  variance,  cr^^,  is  given  in  Table  4.1.  Re-writing  (A.  15)  in  discrete¬ 
time  form  yields  the  following  expression: 

dx,+  i  =  Mrfdx,-  -t  Tc^bc,  -t  (TsJ'a.yfa,  +  -t  cr,^r^,w,4,,  dXo  =  0  (A.  16) 

where 

w^.  is  a  3  X  1  vector  of  independent,  Gaussian  random  variables:  Wa,  ~  N  {0, 1)  and  b^.  is  a 
3x1  vector  of  constant,  random  biases:  bA,.^,  =  b^^,  b^^  ~  N  (o,  cr^^l)- 

The  accelerometer  and  gyro  biases  may  be  augmented  into  the  navigation  state  vector 
as  follows: 

Sxi 

SXi=  bA. 

bo. 

such  that  in  the  absence  of  alignment  errors,  the  discrete-time  navigation  state  error 
dynamics  equation  in  (A.  16)  becomes 

0,.,  0  0 

dx;+i  =  Adx,  +  G  VQw,-,  6x0- U  0,03,9  o-^J  0  (A.  17) 

O3X9  0  ffc  I 

V  L  ^  *0  \) 
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where 
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Appendix  B:  Calibration  for  Bias  and  Drift  Uncertainty 


In  the  conventional  INS  and  in  the  CAI-based  Inertial  Measurement  Unit  (IMU),  the 
accelerometer  and  gyroscope  measurement  errors  at  discrete  time  i  =  1,2,...,  caused 
by  drift  are  white,  Gaussian  noises  Wa,  ~  AfCO,  1)  and  Wg,  ~  AfCO,  1),  respectively 
[15].  Similarly,  the  accelerometer  and  gyroscopes’  unknown,  but  constant,  bias  errors  are 
bA  ~  N(0,crl/)  and  bo  ~  N  (O,  (tIJ)-  The  biases  Ba  and  be  are  3  x  1  vectors  and  their 
uncertainties  are  quantified  by  the  accelerometer  and  gyroscope  specifications,  and 
respectively. 

We  determine  cr„^,  cr^^,  and  by  calibration,  so  that  the  conventional,  navigation 
grade  INS  yields  a  positional  error  of  1  (km/  Vhr)  and  the  CAI-based  IMU  in  theory  yields 
a  position  error  of  5  (m/  Vhr).  This  entails  a  covariance  analysis. 

The  covariance  of  the  navigation  state’s  error  satisfies  the  Lyapunov  equation  as 
follows  [13]: 

=  APf^V  +  GQG^,  i  =  0, 1, ...,  N  -  1  (B.l) 


where  Q  is  the  accelerometers’  and  gyroscopes’  drift  uncertainties  matrix 

0 


Q 


0 


and  N  is  such  that  one  hour  is  considered.  Assuming  no  alignment  errors,  the  covariance 
matrix  is  initialized  as  follows: 


^9x9 

^9x3 

O9X3 

p(^^)  _ 

■^0  “ 

^3x9 

03X3 

^3x9 

03X3 

and  it  is  required  that  the  x  position  uncertainty  of  the  conventional  INS  after  one  hour  be 
specified  as  follows: 

V(P™  ). ,  =  1000  m 
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Note  that  there  is  a  linear  relationship  between  the  sensors’  drift  induced  measurement 
error  variances,  and  crlj  the  variance  and  of  sensor  biases;  and  the  uncertainty 
in  the  aircraft  final  x  position: 

,  =  acrl  +  fieri  +  a.crl  +  fi.crl^ 

The  coefficients  a,  fi,  a^,  and  fi^  are  determined  by  sequentially  setting  all  but  one 
uncertainty  to  zero  in  a  sequential  manner  and  solving  equation  (B.l)  four  times  as  follows: 

O'  —  (Pseoo  w)^  I  ’  ^Vd  ~  ^dd  ~  O’  ^bA  ~  0,  Cfhg  —  0 

P  =  ’  O-v,  =  0,  (Ted  =  1’  =  0’  Cr*G  =  0 

(^h  —  (Pseoow)^  I  5  ^Vd  ~  0,  (TOd  —  0’  ~  ^bc  —  0 

Pb  =  (Pse^y , ,  ’  O-v,  =  0,  (TOd  =  0,  (TbA  =  0,  CT^e  =  1 

The  uncertainty  standard  deviations  are  then  determined  from  the  following  equations: 

’  (^Od  = 

1  CTbc  — 

where  (P^^qq  v)i  i  required  ;c  position  variance  after  a  one  hour  flight. 

For  the  conventional  free  INS,  where  6t  =  0.1  seconds  (N  =  10),  we  require 
(^3«30  v)i  1  “  ^0^  (m^),  such  that  after  one  hour  the  x-axis  position  uncertainty  standard 
deviation  is  1  (km).  The  calculated  measurement  variances,  crl  and  cr^^,  which  account  for 
the  acceleration  and  gyro  drifts  and  the  residual  bias  uncertainties  crl  and  cr^^,  are  given 
in  Table  4.1  for  both  the  INS  and  CAI  IMUs. 

The  CAI  IMU’s  discrete-time  navigation  state  error  covariance  equations  are  the  same 
as  (B.l),  with  the  required  (P36oov)i  j  theoretical 

position  uncertainty  standard  deviation  is  5  (m). 
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